clear
rng default

%% Load useful matrices
load PYX_cond 
load U_CS
x=4;
PYX_cond_x=PYX_cond(x,:);

%% Useful parameters 
n=8^5;
p=n*(n-1)/2;
m=8*10^6;
m_restricted_all=[10^5 m]; 
random_indices_pairs_all=cell(size(m_restricted_all,2),1);
random_indices_pairs_all{end} = randperm(p, m).';
for t=1:size(m_restricted_all,2)-1
    random_indices_pairs_all{t}=randperm(p, m_restricted_all(t)).';
end

Ueq=0;
Uineq=0;
tol=10^(-4);

param_MOSEK.MSK_IPAR_LOG = 0; 
param_MOSEK.MSK_IPAR_INTPNT_BASIS = 'MSK_BI_NEVER';
param_MOSEK.MSK_IPAR_NUM_THREADS = 1;

%% Construct grid
u0grid=0;
u1grid=U_CS(4,1);
gran=60; 
% try different values of gran and different sizes of grid hypercube as discussed in Appendix C of the paper (need powerful cluster to run the code for higher values of gran and larger sizes of hypercube)
u2grid=linspace(-15,15,gran); 
u3grid=linspace(-15,15,gran); 
u4grid=linspace(-15,15,gran); 
u5grid=linspace(-15,15,gran); 
[ca, cb, cc,cd, ce,cf] = ndgrid(u0grid, u1grid, u2grid, u3grid, u4grid, u5grid);
u0grid=ca(:);
u1grid=cb(:);
u2grid=cc(:);
u3grid=cd(:);
u4grid=ce(:);
u5grid=cf(:);
%% Workers
workers=600; %number of parallel workers
jobs=round(size(u2grid,1)/workers); %number of jobs per parallel worker

